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We present a twofold contribution to the numerical simulation of Lindblad equations. First, 
an adaptive numerical approach to approximate Lindblad equations using low-rank dynamics is 
described: a deterministic low-rank approximation of the density operator is computed, and its 
rank is adjusted dynamically, using an on-the-fly estimator of the error committed when reducing 
the dimension. On the other hand, when the intrinsic dimension of the Lindblad equation is too high 
to allow for such a deterministic approximation, we combine classical ensemble averages of quantum 
Monte Carlo trajectories and a denoising technique. Specifically, a variance reduction method based 
upon the consideration of a low-rank dynamics as a control variate is developed. Numerical tests for 
quantum collapse and revivals show the efficiency of each approach, along with the complementarity 
of the two approaches. 
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INTRODUCTION 

Lindblad equations are notoriously challenging to sim¬ 
ulate numerically. The two categories of approaches are 
deterministic approaches, on the one hand, and Monte- 
Carlo approaches [1], on the other hand. Both categories 
have their pros and cons. In the former category, the sim¬ 
ulation is extremely effective when possible, but a major 
difficulty lies in the high dimensionality of the ambient 
space which drastically limits the applicability. In the 
latter category, dimensionality is not an issue, but the 
intrinsic noise of stochastic simulations affects the qual¬ 
ity of the numerical results. The twofold purpose of this 
article is to present recent advances in either of the two 
categories of approaches. 

In a previous work [2], the first two authors have pre¬ 
sented a possible deterministic approach for the simula¬ 
tion of the Lindblad equation, actually borrowed from 
similar ideas introduced in [3, 4] in the context of quan¬ 
tum filtering. The approach consists in approximating 
the evolution of the n x n density matrix p solution to 
the differential Lindblad equation using a reduced dy¬ 
namics on the set of density matrices of some fixed rank 
TO <C n. This reduced dynamics is obtained by taking 
the orthogonal projection of onto the tangent space 
to this set of rank-m matrices . The clear limitation of 
the approach lies in the fact that many practical prob¬ 
lems are not reducible to a low-rank approximation, and 
further that, even when it is the case, the intrinsic di¬ 
mensionality of the reduced dynamics is not necessarily 
known beforehand and may vary in time. So the ques¬ 
tion of adjusting on-the-fly the dimensionality to of the 
low rank dynamics immediately arises. As our first con¬ 
tribution in the present article, we describe below an 
adaptive low-rank simulation, the purpose of which is 


to significantly extend the applicability of the approach 
introduced earlier in [2]. The questions we examine in 
this work enjoy some similarity with questions arising in 
computational quantum chemistry, typically for multi¬ 
configuration time-dependent Hartree and Hartree-Fock 
equations [5, 6]. 

For problems definitely not amenable to determinis¬ 
tic simulation because of their prohibitively high dimen¬ 
sionality (which is indeed the case for many practically 
relevant problems), stochastic approaches are in order, 
see [1, 7-9]. Although the low-rank dynamics no longer 
adequately represents the system, a reduced model, sim¬ 
ulated deterministically, can however serve as a useful 
tool for the stochastic simulation of the high dimensional 
system. We employ the reduced dynamics as a control 
variate within a variance reduction method applied to 
the full, high dimensional stochastic system. Our second 
contribution is to demonstrate the efficiency of such a 
variance reduction method. 


ADAPTIVE LOW-RANK APPROXIMATION 

We consider throughout this article a Lindblad equa¬ 
tion with, for simplicity (and this is by no means a limi¬ 
tation of our methods), a single decoherence operator L, 

^p = -i[H,p\- Lp + pL"' L) + LpL\ (1) 

where p is a n x n non-negative Hermitian matrix with 
Tr (p) = 1, H is & n X n Hermitian matrix and L is 
a n X n matrix. The reduced dynamics derived in [2] 
approximates, for n large, the above dynamics on the set 
of non-negative Hermitian matrices of rank to, to being 
an integer presumably much smaller than n. To make the 
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approximation explicit, one introduces a system of two 
coupled differential equations for U and a corresponding 
to the generic decomposition Plr = U(jU\ where a is 
a, m X m strictly positive Hermitian matrix, U a n x 
m matrix with WU = 1^, and Im denotes the m x m 
identity matrix. That system reads as: 

= -iHU 
dt 

+ (I„ - UU"^) {-^L^LU + LUaU^L'^Ua-'^) , 

( 2 ) 

^cr = -^{U^L^LUa + aU^L^LU) + U'^LUaU'^L^U 
+ iTr(Lt(I„-C/[/t)LC/aC/t)l„. (3) 

Notice that H only appears in (2) and not in (3), a fact 
that is particularly appropriate when H dominates L, in 
which case (3) may be understood as a slow evolution 
as compared to the dynamics (2). Then the projection 
of the original Lindblad dynamics (1) onto the tangent 
space to the set Vm of density matrices of rank m takes 
the explicit form 


where Q = , with V € M", is the one-dimensional 

projector associated to the one dimension added. Denot¬ 
ing by 

G = (In — Pplr)-^Plr.£'^(I„ — Pplr)) 

a straightforward calculation yields 

Tr^ (£j^+Q(p,n)) = ||(Ir - Q)G(ln - Q)f 

+ ;^^Tr2 ((In - Q)G). 

One may then prove (see the details in [10]) that the 
directions 

V = Argminyj_[;_||y||^;^||£p^^^_|_Q(pLR)|| (5) 

minimizing the projection error are the eigenvectors of 
the symmetric matrix G associated to its largest eigen¬ 
value. The matrix G being of large size n, determining 
its largest eigenvector is challenging computationally. To 
this end, we notice that 

V^GV=\\V\lr.-Pp,jLUy^r, 


^^Plr — a PlrA^P) T PPlrP^ 

~ (In ~ PpLR)ApLRA^(In ~ ^Plr) 

I Tr(^PLR.^^(Iti —f*PLR)) p (A\ 

where the orthogonal projection on the image of Plr, 
PpLR = UW, only depends on Plr- Notice that the right¬ 
most term of (4) allows Tr (plr) to be preserved in time. 
In the sequel, we denote by C the right-hand side of (1), 
that of (4), and by = £ — £ll. 

We have described in details in [2] how system (2)- 
(3) may be efficiently simulated and then provides an 
accurate approximation of (1) in the case when the rank 
can be actually reduced. In that work, the rank m was 
prescribed beforehand. Our purpose here is to explain 
how the approach can be amended so as to allow for 
a dynamical adaptation of the rank m of the reduced 
system. 

The adaptation we suggest is based on the evaluation, 
and update, of the projection error 

£ (Plr) = (In ~ ^pLR)ApLRA^(In ~ PpBu) 

Tr (LpRRLt(I„ - ^ 

m 

committed when replacing (1) by (4). 

This error may be reduced upon adding one dimension 
to the m-dimensional subspace Im{U) associated to the 
projector Pp^^. The best possible such dimension to add 
is that for which the projection error is minimal. The 
rank-m projector Pp^^ = UW G is modihed into 

the rank-m-l-1 projector £pLR+Q = UW G ]jmxm 


and that the range of (I„ — Pp^^)LU is of dimension less 
or equal to m and is orthogonal to the range of Plr- Thus, 
denoting by {'pj)o<j<r an orthonormal basis of the range 
of (I„ — Pp^^)LU with r < m, it is sufficient to consider 
V as linear combination of the (fj to get an eigenvector 
of G with largest eigenvalue: V = where the r- 

dimensional vector v of component Vj corresponds to the 
eigenvector with largest eigenvalue of 

K = <l>t(I„ - 

with $ the n X r matrix formed by the r vectors (f>j. 
Since K is of size r < m, this provides an effective manner 
to determine the optimum in (5). 

We have performed a comprehensive series of test of 
the approach. The practical implementation of the dy¬ 
namical adaptivity of the rank is performed as follows. 

We denote by 0{pbr) = j fix a maximal an- 

II^"(plr)II 

gular error O^ax and update the rank 

• increasing m by 1, when ^(plr) > 0max and then 
complementing U via the solution of (5); 

• reducing to by 1, when the smallest eigenvalue A^m 
of (T is such that d(pLj^) -t- \min P 

As an illustrative example, we consider a qubit resonantly 
coupled to a quantized harmonic damped oscillator: 

^P= ^[a^cr.-acr4.,p]-«:(np/2-kpn/2-apa^), (6) 

with Do > 0 the vacuum Rabi pulsation, a the pho¬ 
ton annihilation operator, cr. the qubit lowering operator 
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(cr. = |g)(e|), (T_|_ = aJ, n = a^a the photon-number op¬ 
erator and 1 /k > 0 the oscillator damping time. In Fig¬ 
ure 3.20, page 156 of [9], numerical simulations of quan¬ 
tum collapse and revivals are presented, for /c = 0, when 
the qubit is initially in the excited state |e) and the os¬ 
cillator in a coherent state with h = 15 photons. We 
consider here the same system with k = llo/hOO, a value 
small compared to Hq. This corresponds to a photon 
life-time I/k around 10 times as large as the revival time 
Tr = In this simple case, we can perform the 

full-rank simulation with high precision. This simulation 
illustrated on Figure (1) provides us with a reference cal¬ 
culation. As shown in [2], a constant rank of 4 is suf¬ 
ficient to compute accurately the solution t ^ pt for t 
between 0 and 2Tr. For intermediate values of t, larger 
than 2Tr but not excessively larger, the quantum state p 
gets more and more mixed. For t very large, p becomes 
again pure, since its limit for t i—>■ -l-oo is the lowest en¬ 
ergy state (qubit in the ground state \g), oscillator with 
zero photon). This behavior, illustrated by the numerical 
simulations of Figure 2, is well captured by our adaptive 
approach. 


DENOISED MONTE-CARLO APPROACH 




In the case when the Lindblad equation (1) is gen¬ 
uinely high-dimensional, the model reduction previously 
described is likely to be either ineffective or inaccurate, 
while the direct integration of the equation is out of 
reach. In this situation, the classical approach is to use 
a Monte-Carlo sampling. One derives a stochastic dy¬ 
namics on the wave function [ipt) such that the density 
p{t) = M {\tjjt) {tjjtl) constructed from Itljt) (E stands for 
expectation value, i.e. ensemble average) solves the Lind¬ 
blad equation (1). In practice, M independent trajecto- 



with A: = 1,..., M, are then simulated using 


that stochastic dynamics and the estimator of the mean 


M 


PMC (A) = 




(k) 




( 7 ) 


is used as an approximation of p{t). The practical diffi¬ 
culty of Monte-Carlo approaches is, as briefly mentioned, 
above, the variance, that is, the noise intrinsically present 
in the approach. 

In principle, there are infinitely many dynamics on \ipt) 
that are consistent with the Lindblad dynamics. Interest¬ 
ingly, a straightforward calculation from (7) shows that 

E (Tr ((pMC - pf) ) = — (8) 

The variance of the estimator of the mean pMC is 
therefore independent of the specihc unravelling choice, 
namely the stochastic dynamics set on |'(/it), provided that 


FIG. 1. (color online) High precision full-rank numerical sim¬ 
ulation of f I— p{t) for the quantum collapse and revivals of a 
qubit resonantly coupled to a slightly damped quantum har¬ 
monic oscillator governed by (6). Top plot: the solid curve 
corresponds to the excited population (e|p|e) versus time; the 
normalized time corresponds to t/Tr where Tr is the revival 
time; the oscillator damping time is around lOTr; the initial 
state po corresponds to the qubit in excited state |e) and os¬ 
cillator in a coherent state with n = 15 photons; a truncation 
of the number of photons to a maximum of 2h = 30 yields 
an underlying Hilbert space of dimension 62. Bottom plot: 
zoom of top plot for a normalized time between 0 and 10. 


state remains normalized. This is easily seen in the 
proof of (8) . This property is, of course, a remarkable pe¬ 
culiarity of the present context. And the discretization 
in time of the process can of course slightly affect that 
property. The classical unravelling choice (see [11-13]) is 
the Wiener process defined by 

d\%l^t) = Di{\%l;t)) dt + D2{\i!t))dWt., (9) 

with the drift term 

-Di(l^)) = |'0) + T>°(|i/))) where 

D%m = \{{L + L^)wL- L^L - i(L + L%^'^ \^) , 

( 10 ) 

and the diffusion 

D2m)=(^L-^-{L + L%^y^}. ( 11 ) 
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normalized time 



FIG. 2. (color online) Comparison between the full-rank tra¬ 
jectory 1 1 —>■ pit) illustrated on figure 1 with the adaptive rank 
trajectory t i—>■ pLB.{t) governed by (2-3). Top plot: the solid 
curve corresponds to -^/Tr ((p — Plr)^) versus time with scale 
on the left. The dashed curve corresponds to the adaptive 
rank m with scale on the right; the initial low rank state Plr 
coincides with po; the initial rank m is thus set to one; it 
evolves according to the maximal angular error Omax = . 

Bottom plot: zoom of top plot for a normalized time between 
0 and 10. 

We have used the notation (A)|^) = (^/>|^|'!/;). We empha¬ 
size that other choices of dynamics on \tp), all consistent 
with the Lindblad equation through p{t) = E(|V^t) 
could be made. In particular, Poisson processes could 
be considered instead of Wiener processes (a choice that 
can be seen as more natural given the applications ad¬ 
dressed in the present work where photons are emitted or 
absorbed). In any event, given the property (8) and as¬ 
suming that the Poisson process remains normalized, the 
variance remains identical. We indeed double-checked 
that, in actuality, using Poisson processes does not bring 
further practical variance reduction, see [10] for more de¬ 
tails. 

Variance is an issue for the numerical simulation and 
affects the accuracy of the results. It is thus desirable 
to come up with further, dedicated variance reduction 
approaches that may reduce the computational cost at 
accuracy fixed, or improve the accuracy for a given com¬ 


putational cost. An approach that has proved effective 
in many engineering sciences for reducing variance of 
Monte-Carlo simulations is that of control variate, see 
e.g. [14, page 54]. In short, the approach consists in con¬ 
currently simulating the original system under consid¬ 
eration and a system correlated to that original system 
so as to minimize the variance in the simulation of the 
former system. Intuitively, the approach works by ’’can¬ 
cellation” of the noise because the same random draws 
are used for both systems. More specifically, we consider 
here as control variate the low-rank approximation Plr 
of previous section. Even though that low-rank system is 
not a correct approximation of the original system (which 
we have deliberately assumed here high-dimensional), it 
is sufficiently correlated to that system to provide an ef¬ 
ficient variance reduction. Practically, we construct an 
(’’Control-Variate”) estimated density depending on the 
adjustable scalar parameter A, 

Pcv = pMC + -^(plr — Pmclr)i (12) 

as a combination of the original Monte-Carlo estimated 
density pmc (constructed from the simulation of (7)- 
(9)), the estimated density pmclr constructed from a 
low-rank dynamics (” Monte-Carlo-Low-Rank”, see be¬ 
low (13)), and the density plr obtained upon solving 
the corresponding low-rank Lindblad equation. Since by 
construction 

E[pMCLR(i)] = Plr(1), 

the approximation method (12) is unbiased (taking the 
expectation of both sides of (12) yields p = E[pcv] = 
E[pmc])- The scalar parameter A is adjusted so as to 
minimize the variance E ^Tr ^(pcv ~ p)^^ ^ i a polyno¬ 
mial in A of degree 2, in such a way that 

E (Tr ((pcv - P)')) « E (Xr ((pMC - p)')) 

so that the simulation of pcv is eventually more effec¬ 
tive than that of pmc- The more correlated the reduced 
model and the original model, the closer A to one, the 
smaller this variance and thus the more efficient the de- 
noising. In passing, we notice that, although this will not 
be the case in the actual numerical experiments we per¬ 
form, the low-rank dynamics plr( 1) could itself be chosen 
with an adaptive rank, as in the previous section. Like¬ 
wise, we could pick as control variate another dynamics 
than the low-rank dynamics, if a more convenient one is 
available. 

The remaining question is to derive a stochas¬ 
tic dynamics on a wave function J'^lr) such that 
EdV'LR) (i/’lrI) = Plr) given that the dynamics of Plr is 
known since easy to compute via [U, a) solutions of (2) 
and (3). For this purpose, it is a natural idea to seek Ji/'lr) 
under the form JV'lr) = U \v) where the reduced wave 
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function \v) € K™ is a stochastic process to be deter¬ 
mined. It can be shown (see [10]) that the correct dy¬ 
namics to consider reads 


d IV’lr) (IV’RR))dt + 

+ (I„ - Pp^J + LUaU^L^Ua-^U^^ dt 


Tr((I„ 


-PpLR)^PLRi^) 

2m 


Ua IV'lr) dt, 


(13) 


with Di and D 2 defined in (10) and (11). 


Since pmclr = jj J2k=i 


^£>) {i’S 


a simple com¬ 


putation, exploiting the fact that both and 

are normalized , shows that 




M E (Xr ((pcv - P)")) = (1 - Tr (p^)) 

+ 2 (e (K^I^rr)!") - Tr(ppR,)) A + 1 - Tr (p^) . 
Thus the optimal denoising choice for A reads 


A = 


Tr( pPlr) - E 


l-Tr(pR^2) 


(14) 


In practice, and also in the simulations of Figure (3), the 
adjustable parameter A is given by (14) where the p is 
replaced by pMC (we do not have access to p itself) and 

where E ^I(V^IV^lr) is replaced by the estimated mean 

Figure (3) illustrates the interest of this denoising 
method for the Lindblad system simulated in Figure (2). 
We take M = 400 trajectories. The dynamics of I'^lr) is 
based on a low-rank approximation Prr of constant rank 
m — 2. On the bottom plot of Figure 2, we observe 
that, at normalized time 2, an accurate approximation 
of p must be of rank 10 or larger. Nevertheless, Fig¬ 
ure 3 indicates that, at normalized time 2, a reduction 
of the standard deviation around 50% is obtained with 
pcv instead of pmc- Such a variance reduction corre¬ 
sponds, with a classical quantum Monte-Carlo method, 
to an increase of the number of trajectories by a factor 
4. Such a gain confirms the definite interest of combin¬ 
ing deterministic low-rank approximations with quantum 
Monte-Carlo trajectories for the numerical simulation of 
high-dimensional Lindblad equation. 
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normalized time 

FIG. 3. (color online) Variance reduction corresponding to the 
quantum collapse and revivals in the bottom plot of Figure 2. 
The solid black curve corresponds to -y/Tr ((p — pcv)^), the 
dashed blue one to -y/Tr ((p — Pmc)^) and the dotted red one 
to A. The number of trajectories M is set to 400. The low- 
rank approximation Plr used for the control variate governed 
by (13) is maintained at the constant rank m = 2. 
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